# Manuscript: Signaling Resolve Through Credit-Claiming
# Author: Ilayda B. Onder, The Pennsylvania State University, Department of Political Science
# Replication File for International Interactions
# April 8, 2023

# This analysis is run using R version 4.2.1 (2022-06-23) for Mac
# The codes in this replication file produces Table 1, Table 2, Table 3 and Figure 1
# The sections in this replication file are in order. Please run them one by one. 

# The packages needed to be installed to run the analysis:
install.packages("stargazer") # for publication quality tables
install.packages("sampleSelection") # for Heckman selection models
install.packages("pscl") # for Zero-inflated negative binomial models
install.packages("ggplot2") # for graphs
install.packages("dplyr") # for dataframe manipulation

# Load packages
library(stargazer)
library(sampleSelection)
library(pscl)
library(ggplot2)
library(dplyr)

# WARNING: Please set the working directory to the current Source File Location by clicking Session > Set Working Directory > To Source File Location

# Load analysis datasets
part1 = get(load("Analysis Data Part 1.RData")) # This dataset is used to run the models presented in Table 1
part2 = get(load("Analysis Data Part 2.RData")) # This dataset is used to run the models presented in Tables 2 and 3, and Figures 1, 2, and 3
rm(dat, dt) # remove duplicate dataframes

# Table 1: Logit Models of Retaliatory and Conciliatory Government Response against Militant Groups, 1998-2012----
# Use the dataframe "part1" for replicating Table 1
## DV: Retaliation
### Model 1
fit1 <- glm(stick ~ 
              claimed_log_lag1 + # claimed attacks
              unclaimed_log_lag1 + # unclaimed attacks
              ucdpbd_log + # battle deaths
              nonterr_casualties_log + # non-militant casualties
              size_rec + # group size
              ethn + # ethno-nationalist group
              reli + # religious group
              left.x, # leftist group
            family = binomial, data = part1)
### Model 2
fit2 <- glm(stick ~ 
              claimed_log_lag1 + # claimed attacks
              unclaimed_log_lag1 + # unclaimed attacks
              ucdpbd_log + # battle deaths
              nonterr_casualties_log + # non-militant casualties
              size_rec + # group size
              ethn + # ethno-nationalist group
              reli + # religious group
              left.x + # leftist group
              state_sponsor + # state sponsor
              terrctrl + # territorial control
              socsvcs + # social service provision
              holdoffice, # legal political office
            family = binomial, data = part1)
### Model 3
fit3 <- glm(stick ~ 
              claimed_log_lag1 + # claimed attacks
              unclaimed_log_lag1 + # unclaimed attacks
              ucdpbd_log + # battle deaths
              nonterr_casualties_log + # non-militant casualties
              size_rec + # group size
              ethn + # ethno-nationalist group
              reli + # religious group
              left.x + # leftist group
              state_sponsor + # state sponsor
              terrctrl + # territorial control
              socsvcs + # social service provision
              holdoffice + # legal political office
              gd_ptsa_inf + # state human rights violations
              fh_ipolity2_inf + # state regime type
              gdp_per_capita_un_log + # state gdp per capita logged
              logwdi_pop, # state population logged
            family = binomial, data = part1)
## DV: Conciliation
### Model 4
fit4 <- glm(carrot ~ 
              claimed_log_lag1 + # claimed attacks
              unclaimed_log_lag1 + # unclaimed attacks
              ucdpbd_log + # battle deaths
              nonterr_casualties_log + # non-militant casualties
              size_rec + # group size
              ethn + # ethno-nationalist group
              reli + # religious group
              left.x, # leftist group
            family = binomial, data = part1)
### Model 5
fit5 <- glm(carrot ~ 
              claimed_log_lag1 + # claimed attacks
              unclaimed_log_lag1 + # unclaimed attacks
              ucdpbd_log + # battle deaths
              nonterr_casualties_log + # non-militant casualties
              size_rec + # group size
              ethn + # ethno-nationalist group
              reli + # religious group
              left.x + # leftist group
              state_sponsor + # state sponsor
              terrctrl + # territorial control
              socsvcs + # social service provision
              holdoffice, # legal political office
            family = binomial, data = part1)
### Model 6
fit6 <- glm(carrot ~ 
              claimed_log_lag1 + # claimed attacks
              unclaimed_log_lag1 + # unclaimed attacks
              ucdpbd_log + # battle deaths
              nonterr_casualties_log + # non-militant casualties
              size_rec + # group size
              ethn + # ethno-nationalist group
              reli + # religious group
              left.x + # leftist group
              state_sponsor + # state sponsor
              terrctrl + # territorial control
              socsvcs + # social service provision
              holdoffice + # legal political office
              gd_ptsa_inf + # state human rights violations
              fh_ipolity2_inf + # state regime type
              gdp_per_capita_un_log + # state gdp per capita logged
              logwdi_pop, # state population logged
            family = binomial, data = part1)
## Produce Table 1
stargazer(fit1, fit2, fit3, fit4, fit5, fit6, 
          type = "text")



# Table 2: Heckman Selection Models of Militant Credit-Claiming, 1998-2016----
# Use the dataframe "part2" to replicate Table 2
part2$claimedratio = part2$claimed / part2$all # create the ratio of claimed attacks to the total number of attacks
part2$stage1 = ifelse(part2$all > 0, 1, 0) # 1st stage is whether there was at least one attack, create stage1 variable using the number of all attacks
part2$stage2 = ifelse(part2$claimedratio == 1, 1, 0) # 2nd stage is whether all attacks were claimed, create stage1 variable using the ratio of claimed attacks
part2$rivaldum = ifelse(part2$num_rivals == 0, 0, 1) # create a dummy variable for the presence of rival groups using the number of rivals variable

# Model 7
s1 <- selection(stage1 ~ mul_bases +
                  duration,
                stage2 ~ 
                  duration + 
                  all_log +
                  nonterr_casualties_log + diversity + shr_trans + rel + 
                  splinterhist + civi,
                data = part2)
# Model 8
s2 <- selection(stage1 ~ mul_bases +
                  duration,
                stage2 ~ 
                  duration*splinterhist + 
                  all_log +
                  nonterr_casualties_log + diversity + shr_trans + rel,
                data = part2)
# Model 9
s3 <- selection(stage1 ~ mul_bases +
                  duration,
                stage2 ~ 
                  duration*civi +
                  all_log +
                  nonterr_casualties_log + diversity + shr_trans + rel,
                data = part2)
# Model 10
s4 <- selection(stage1 ~ mul_bases + 
                  fate_leader_dummy,
                stage2 ~ fate_leader_dummy + fate_leader_dummy_lag1 + fate_leader_dummy_lag2 + fate_leader_dummy_lag3 + 
                  splinterhist + civi +
                  all_log + 
                  nonterr_casualties_log + diversity + shr_trans + rel,
                data = part2)
# Model 11
s5 <- selection(stage1 ~ mul_bases +
                  duration,
                stage2 ~ 
                  duration + 
                  all_log +
                  nonterr_casualties_log + diversity + shr_trans + rel + 
                  splinterhist + civi +
                  lead_hierarch + rivaldum,
                data = part2)
# Model 12
s6 <- selection(stage1 ~ mul_bases +
                  duration,
                stage2 ~ 
                  duration*splinterhist + 
                  all_log +
                  nonterr_casualties_log + diversity + shr_trans + rel +
                  lead_hierarch + rivaldum,
                data = part2)
# Model 13
s7 <- selection(stage1 ~ mul_bases +
                  duration,
                stage2 ~ 
                  duration*civi +
                  all_log +
                  nonterr_casualties_log + diversity + shr_trans + rel +
                  lead_hierarch + rivaldum,
                data = part2)
# Model 14
s8 <- selection(stage1 ~ mul_bases + 
                  fate_leader_dummy,
                stage2 ~ fate_leader_dummy + fate_leader_dummy_lag1 + fate_leader_dummy_lag2 + fate_leader_dummy_lag3 + 
                  splinterhist + civi +
                  all_log + 
                  nonterr_casualties_log + diversity + shr_trans + rel +
                  lead_hierarch + rivaldum,
                data = part2)
# Produce Table 2
stargazer(s1, s2, s3, s4, s5, s6, s7, s8, 
          type = "text", 
          order = c("duration", "splinterhist", "civi1", "duration*splinterhist", "duration:civi1", # order the variables 
                    "fate_leader_dummyDecapitation", "fate_leader_dummy_lag1Decapitation", "fate_leader_dummy_lag2Decapitation", "fate_leader_dummy_lag3Decapitation",
                    "all_log", "nonterr_casualties_log", "diversity", "shr_trans", "rel1", "lead_hierarch1", "rivaldum"),
          covariate.labels = c("Group duration (in years)", # fix variable names
                               "Duration*Splinter group",
                               "Duration*Civilian targets",
                               "Splinter group",
                               "Civilian targets",
                               "Decapitation",
                               "Decapitation (L1)",
                               "Decapitation (L2)",
                               "Decapitation (L3)",
                               "Total number of attacks",
                               "Number of non-militant casualties",
                               "Attack diversity",
                               "Share of transnational attacks",
                               "Religious fundamentalist group",
                               "Hierarchical leadership",
                               "Rival groups",
                               "Constant"),
          column.labels = c("Duration", "Splinter", "Target", "Decapitation", "Duration", "Splinter", "Target", "Decapitation"), # create labels for each model
          model.numbers = T) # include model numbers




# Table 3: Zero-Inflated Negative Binomial Models of Militant Credit-Claiming, 1998-2016----
# Use the dataframe "part2" to replicate Table 3
# Model 15
z1 <- zeroinfl(claimed ~ 
                 duration + 
                 all_log +
                 nonterr_casualties_log + diversity + shr_trans + rel + 
                 splinterhist + civi |
                 duration,
               data = part2, 
               dist = "negbin") # specify distribution of negative binomial
# Model 16
z2 <- zeroinfl(claimed ~ 
                 duration*splinterhist + 
                 all_log +
                 nonterr_casualties_log + diversity + shr_trans + rel |
                 duration,
               data = part2, 
               dist = "negbin") # specify distribution of negative binomial
# Model 17
z3 <- zeroinfl(claimed ~ 
                 duration*civi +
                 all_log +
                 nonterr_casualties_log + diversity + shr_trans + rel |
                 duration,
               data = part2, 
               dist = "negbin") # specify distribution of negative binomial
# Model 18
z4 <- zeroinfl(claimed ~ 
                 fate_leader_dummy + fate_leader_dummy_lag1 + fate_leader_dummy_lag2 + fate_leader_dummy_lag3 + 
                 splinterhist + civi +
                 all_log + 
                 nonterr_casualties_log + diversity + shr_trans + rel |
                 fate_leader_dummy,
               data = part2, 
               dist = "negbin") # specify distribution of negative binomial
# Model 19
z5 <- zeroinfl(claimed ~ 
                 duration + 
                 all_log +
                 nonterr_casualties_log + diversity + shr_trans + rel + 
                 splinterhist + civi +
                 lead_hierarch + rivaldum |
                 duration,
               data = part2,
               dist = "negbin") # specify distribution of negative binomial
# Model 20
z6 <- zeroinfl(claimed ~ 
                 duration*splinterhist + 
                 all_log +
                 nonterr_casualties_log + diversity + shr_trans + rel +
                 lead_hierarch + rivaldum |
                 duration,
               data = part2, 
               dist = "negbin") # specify distribution of negative binomial
# Model 21
z7 <- zeroinfl(claimed ~ 
                 duration*civi +
                 all_log +
                 nonterr_casualties_log + diversity + shr_trans + rel +
                 lead_hierarch + rivaldum |
                 duration,
               data = part2, 
               dist = "negbin") # specify distribution of negative binomial
# Model 22
z8 <- zeroinfl(claimed ~ 
                 fate_leader_dummy + fate_leader_dummy_lag1 + fate_leader_dummy_lag2 + fate_leader_dummy_lag3 + 
                 splinterhist + civi +
                 all_log + 
                 nonterr_casualties_log + diversity + shr_trans + rel +
                 lead_hierarch + rivaldum |
                 fate_leader_dummy,
               data = part2, 
               dist = "negbin") # specify distribution of negative binomial
# Produce Table 3
stargazer(z1, z2, z3, z4, z5, z6, z7, z8,
          type = "text", 
          order = c("duration", "splinterhist", "civi1", "duration*splinterhist", "duration:civi1", # reorder variables
                    "fate_leader_dummyDecapitation", "fate_leader_dummy_lag1Decapitation", "fate_leader_dummy_lag2Decapitation", "fate_leader_dummy_lag3Decapitation",
                    "all_log", "nonterr_casualties_log", "diversity", "shr_trans", "rel1", "lead_hierarch1", "rivaldum"),
          covariate.labels = c("Group duration (in years)", # fix variable names
                               "Duration*Splinter group",
                               "Duration*Civilian targets",
                               "Splinter group",
                               "Civilian targets",
                               "Decapitation",
                               "Decapitation (L1)",
                               "Decapitation (L2)",
                               "Decapitation (L3)",
                               "Total number of attacks",
                               "Number of non-militant casualties",
                               "Attack diversity",
                               "Share of transnational attacks",
                               "Religious fundamentalist group",
                               "Hierarchical leadership",
                               "Rival groups",
                               "Constant"),
          column.labels = c("Duration", "Splinter", "Target", "Decapitation", "Duration", "Splinter", "Target", "Decapitation"), # add labels for each model
          model.numbers = T) # include model numbers








# Figure 1: Trends in Credit-Claiming Behavior and Group Survival----
# Use the dataframe "part2" to replicate Figure 1
## Panel A: Annual trends in militant attacks
desc = aggregate(claimed ~ year, part2, sum) # aggregate the number of claimed attacks over year
desc2 = aggregate(unclaimed ~ year, part2, sum) # aggregate the number of unclaimed attacks over year
desc = merge(desc, desc2, by = "year") # merge the number of claimed and unclaimed attacks into the same dataframe

plot(desc$year, desc$claimed, type="o", col = gray(0, .8), pch="o", lty=1, lwd = 3, ylim=c(0,4500), # plot claimed attacks over year
     xlab = "Year", ylab = "Number of attacks")
points(desc$year, desc$unclaimed, col = gray(0, .4), pch="*", cex = 3) # add unclaimed attacks onto the same plot
lines(desc$year, desc$unclaimed, col = gray(0, .4),lty=2, lwd = 3)
legend(2000,4000,legend=c("Claimed attacks", "Unclaimed attacks"), col=c(gray(0, .8), gray(0, .4)), # create the legend
       pch=c("o","*"),lty=c(1,2), lwd = 3, ncol=1)

## Panel B: Survival rates of militant groups
dursplit = part2 %>% # only keep the last observation in the group duration variable, which signifies the maximum group duration for each group
  group_by(gname) %>%
  filter(duration == max(duration))
dursplit = subset(dursplit, !is.na(fate_leader_dummy)) # drop observations whose decapitation variable is NA

hist(dursplit$duration[dursplit$fate_leader_dummy == "No Decapitation"], # Plot the histogram for groups that never experienced decapitation
     main = "",
     xlab = "Group duration (in years)",
     ylab = "Frequency",
     breaks = 10,
     xlim = c(0, 50),
     col = gray(0, .4))
hist(dursplit$duration[dursplit$fate_leader_dummy == "Decapitation"], # add the histogram for groups that experieced decapitation to the same plot
     breaks = 10,
     col = gray(0, .8),
     add = T)
legend(25,200,legend=c("Non-decapitated groups", "Decapitated groups"), col=c(gray(0, .4),gray(0, .8)), # add legend
       lty=c(1,1), lwd = 10, ncol=1)







# END----






